function varargout = BTL_XLDLDPT(varargin)
% BTL_XLDLDPT MATLAB code for BTL_XLDLDPT.fig
%      BTL_XLDLDPT, by itself, creates a new BTL_XLDLDPT or raises the existing
%      singleton*.
%
%      H = BTL_XLDLDPT returns the handle to a new BTL_XLDLDPT or the handle to
%      the existing singleton*.
%
%      BTL_XLDLDPT('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in BTL_XLDLDPT.M with the given input arguments.
%
%      BTL_XLDLDPT('Property','Value',...) creates a new BTL_XLDLDPT or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before BTL_XLDLDPT_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to BTL_XLDLDPT_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help BTL_XLDLDPT

% Last Modified by GUIDE v2.5 13-Apr-2013 16:41:30

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @BTL_XLDLDPT_OpeningFcn, ...
                   'gui_OutputFcn',  @BTL_XLDLDPT_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before BTL_XLDLDPT is made visible.
function BTL_XLDLDPT_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to BTL_XLDLDPT (see VARARGIN)

% Choose default command line output for BTL_XLDLDPT
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes BTL_XLDLDPT wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = BTL_XLDLDPT_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


% --- Executes on button press in imgSelect.
function imgSelect_Callback(hObject, eventdata, handles)
% hObject    handle to imgSelect (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
axes(handles.originalImg);
[filename, pathname] = uigetfile('*.bmp','Select an image file');
img = imread(fullfile(pathname,filename));
set(handles.imgLink,'String',fullfile(pathname,filename));
handles.image = img;
guidata(hObject,handles);
imshow(img);


% --- Executes on button press in imgCompress.
function imgCompress_Callback(hObject, eventdata, handles)
% hObject    handle to imgCompress (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if ~isfield(handles,'image')
    handles.image = imread('cameraman.tif');
end

%-------------------JPEG2000-----------------------%
NL = 3; % number of levels of 2D DWT
SamplingFormat = '4:2:0';
A = handles.image;
%A = imread('cameraman.tif');

%change image size to make sure that the image size is divisible by NL
[X,Y,Z] = size(A);
if mod(X , 2^NL)~=0
    Height = floor(X/(2^NL))*(2^NL);
else
    Height = X;
end

if mod(Y,2^NL)~=0
    Width = floor(Y/(2^NL))*(2^NL);
else
    Width = Y;
end

Depth = Z;
clear X Y Z

%level offset -128 (DC level shift
%convert color sytem RGB to YCrCb
% In JPEG2000 jargon, y0 corresponds to the Luma Y,
% y1 to Cb and y2 to Cr.
if Depth == 1
    %gray image
    y0 = double(A(1:Height,1:Width)) - 128; % DC level shift
else
    %color image
    
    A1 = double(A(1:Height,1:Width,:))-128; % DC level shift
    %change color system to YCbCr
    
    y0 = 0.299*A1(:,:,1) + 0.587*A1(:,:,2) + 0.144*A1(:,:,3);
    y1 = -0.16875*A1(:,:,1) - 0.33126*A1(:,:,2) + 0.5*A1(:,:,3);
    y2 = 0.5*A1(:,:,1) - 0.41869*A1(:,:,2) - 0.08131*A1(:,:,3);
    clear A1
    
    %Resize Cr and Cb base on sampling format
    switch SamplingFormat
        case '4:2:0'
            y1 = imresize(y1,[Height/2 Width/2],'cubic');
            y2 = imresize(y2,[Height/2 Width/2],'cubic');
        case '4:2:2'
            y1 = imresize(y1,[Height Width/2],'cubic');
            y2 = imresize(y2,[Height Width/2],'cubic');
    end
end

% Daubechies 9/7 Filters
% Note: The Daubechies 9/7 filter bank has 9 taps for its analysis
% Lowpass and synthesis highpass FIR filters;
% 7 taps for its analysis highpass and synthesis lowpass FIR filters.
% However, in order to obtain subbands that are
% exact multiples of 2^-n, n being a positive integer, the
% number of taps in the filters must be even. That is why you
% will notice the filter lengths to be 10 instead of 9 or 7.
% The extra tap values are zero.
%
%LOD = analysis lowpass filter
%LOR = synthesis lowpass filter
%HID = analysis highpass filter
%HIR = synthesis highpass filter
LOD = [0 0.0267487574 -0.0168641184 -0.0782232665 0.2668641184...
0.6029490182 0.2668641184 -0.0782232665 -0.0168641184...
0.0267487574];
HID = [0 0.0912717631 -0.0575435263 -0.5912717631...
1.1150870525 -0.5912717631 -0.0575435263 0.0912717631 0 0];
%
LOR = [0 -0.0912717631 -0.0575435263 0.5912717631...
1.1150870525 0.5912717631 -0.0575435263 -0.0912717631 0 0];
HIR = [0 0.0267487574 0.0168641184 -0.0782232665 -0.2668641184...
0.6029490182 -0.2668641184 -0.0782232665 0.0168641184...
0.0267487574];
%
% "dwtCoef" cell array has NL x 4 cells.
% In each cell at level NL, the 1st component is the LL
% coefficients,
% 2nd component is HL coefficients, 3rd component LH and
% 4th component HH coefficients.
%
% Do an NL-level 2D DWT of the Y component

dwty0 = cell(NL,4);
for n = 1:NL
if n==1
[dwty0{n,1},dwty0{n,2},dwty0{n,3},dwty0{n,4}] =...
dwt2(y0,LOD,HID,'mode','per');
else
[dwty0{n,1},dwty0{n,2},dwty0{n,3},dwty0{n,4}] =...
dwt2(dwty0{n-1,1},LOD,HID,'mode','per');
end
end
% Quantize the Y DWT coefficients
[dwty0,y0Bits] = jpg2KQuantize(dwty0,NL,1);
%
% Do an NL-level 2D IDWT of the Y component
% The cell array idwtimg has NL cells with each cell
% corresponding to
% the reconstructed LL image at that level.
idwty0 = cell(NL,1);
for n = NL:-1:1
if n==3
idwty0{n}= idwt2(dwty0{n,1},dwty0{n,2},dwty0{n,3},...
dwty0{n,4},LOR,HIR,'mode','per');
else
idwty0{n}= idwt2(idwty0{n+1},dwty0{n,2},dwty0{n,3},...
dwty0{n,4},LOR,HIR,'mode','per');
end
end
SNRy0 = 20*log10(std2(y0)/std2(y0-idwty0{1}));
if Depth > 1
% Do an NL-level 2D DWT of the Cb and Cr components
dwty1 = cell(NL,4);
dwty2 = cell(NL,4);
for n = 1:NL
if n==1
[dwty1{n,1},dwty1{n,2},dwty1{n,3},dwty1{n,4}] =...
dwt2(y1,LOD,HID,'mode','per');
[dwty2{n,1},dwty2{n,2},dwty2{n,3},dwty2{n,4}] =...
dwt2(y2,LOD,HID,'mode','per');
else
[dwty1{n,1},dwty1{n,2},dwty1{n,3},dwty1{n,4}] =...
dwt2(dwty1{n-1,1},LOD,HID,'mode','per');
[dwty2{n,1},dwty2{n,2},dwty2{n,3},dwty2{n,4}] =...
dwt2(dwty2{n-1,1},LOD,HID,'mode','per');
end
end
%
% Quantize the Cb and Cr DWT coefficients
[dwty1,y1Bits] = jpg2KQuantize(dwty1,NL,2);
[dwty2,y2Bits] = jpg2KQuantize(dwty2,NL,3);
%
% % Do an NL-level 2D IDWT of the Cb and Cr components
idwty1 = cell(NL,1);
idwty2 = cell(NL,1);
%
for n = NL:-1:1
if n==3
idwty1{n}= idwt2(dwty1{n,1},dwty1{n,2},dwty1{n,3},...
dwty1{n,4},LOR,HIR,'mode','per');
idwty2{n}= idwt2(dwty2{n,1},dwty2{n,2},dwty2{n,3},...
dwty2{n,4},LOR,HIR,'mode','per');
else
idwty1{n}= idwt2(idwty1{n+1},dwty1{n,2},dwty1{n,3},...
dwty1{n,4},LOR,HIR,'mode','per');
idwty2{n}= idwt2(idwty2{n+1},dwty2{n,2},dwty2{n,3},...
dwty2{n,4},LOR,HIR,'mode','per');
end
end
% Compute the SNR of Cb and Cr components due to quantization
SNRy1 = 20*log10(std2(y1)/std2(y1-idwty1{1}));
SNRy2 = 20*log10(std2(y2)/std2(y2-idwty2{1}));
% Upsample Cb & Cr components to original full size
y1 = imresize(idwty1{1},[Height Width],'cubic');
y2 = imresize(idwty2{1},[Height Width],'cubic');
end
% Do inverse ICT & level shift
if Depth == 1
Ahat = idwty0{1}+ 128;

%subplot(1,2,1), imshow(A), title('Originnal Image');
%subplot(1,2,2), imshow(uint8(Ahat)), title('JPEG2000 Compressed Image');
axes(handles.jpeg2000);
imshow(uint8(Ahat));

disp('Thong so cua JPEG2000')
sprintf('Average bit rate = %5.2f bpp\n',y0Bits/(Height*Width))
SNRy0 = 20*log10(std2(y0)/std2(y0-idwty0{1}-128));

str = sprintf('Average bit rate = %5.2f bpp\n\nSNR(Y) = %4.2f dB\n',y0Bits/(Height*Width),SNRy0);
set(handles.jpeg2000Info,'String',str);
else
Ahat = zeros(Height,Width,Depth);
Ahat(:,:,1) = idwty0{1}+ 1.402*y2 + 128;
Ahat(:,:,2) = idwty0{1}- 0.34413*y1 - 0.71414*y2 + 128;
Ahat(:,:,3) = idwty0{1}+ 1.772*y1 + 128;

%subplot(1,2,1), imshow(A), title('Originnal Image');
%subplot(1,2,2), imshow(uint8(Ahat)), title('JPEG2000 Compressed Image');
axes(handles.jpeg2000);
imshow(uint8(Ahat));

disp('Thong so cua JPEG2000')
sprintf('Average bit rate = %5.2f bpp\n',...
(y0Bits+y1Bits+y2Bits)/(Height*Width))
SNRy0 = 20*log10(std2(y0)/std2(y0-idwty0{1}));
str = sprintf('Bit rate = %5.2f bpp\n\nSNR(Y) %4.2fdB\t\nSNR(Cb) = %4.2f dB\t\nSNR(Cr) =%4.2fdB\n\n',(y0Bits+y1Bits+y2Bits)/(Height*Width),SNRy0,SNRy1,SNRy2);
set(handles.jpeg2000Info,'String',str);
end


%-------------------------compression image jpeg------------------------%

%A = imread('birds.ras'); % read an image
%A = imread('lena.bmp');
[Height,Width,Depth] = size(A);
N = 8; % Transform matrix size
% Limit Height & Width to multiples of 8
if mod(Height,N)~=0
Height = floor(Height/N)*N;
end
if mod(Width,N)~=0
Width = floor(Width/N)*N;
end
A1 = A(1:Height,1:Width,:);
clear A
A = A1;
SamplingFormat = '4:2:0';
if Depth == 1
y = double(A);
else
A = double(rgb2ycbcr(A));
y = A(:,:,1);
switch SamplingFormat
case '4:2:0'
Cb = imresize(A(:,:,2),[Height/2 Width/2],'cubic');
Cr = imresize(A(:,:,3),[Height/2 Width/2],'cubic');
case '4:2:2'
Cb = imresize(A(:,:,2),[Height Width/2],'cubic');
Cr = imresize(A(:,:,3),[Height Width/2],'cubic');
end
end
jpgQstepsY = [16 11 10 16 24 40 51 61;...
12 12 14 19 26 58 60 55;...
14 13 16 24 40 57 69 56;...
14 17 22 29 51 87 80 62;...
18 22 37 56 68 109 103 77;...
24 35 55 64 81 104 113 92;...
49 64 78 87 103 121 120 101;...
72 92 95 98 112 100 103 99];
QstepsY = jpgQstepsY;
Qscale =3.2;
Yy = zeros(N,N);
xqY = zeros(Height,Width);
acBitsY = 0;
dcBitsY = 0;
if Depth >1
jpgQstepsC = [17 18 24 47 66 99 99 99;...
18 21 26 66 99 99 99 99;...
24 26 56 99 99 99 99 99;...
47 66 99 99 99 99 99 99;...
99 99 99 99 99 99 99 99;...
99 99 99 99 99 99 99 99;...
99 99 99 99 99 99 99 99;...
99 99 99 99 99 99 99 99];
QstepsC = jpgQstepsC;
YCb = zeros(N,N);
YCr = zeros(N,N);
switch SamplingFormat
case '4:2:0'
xqCb = zeros(Height/2,Width/2);
xqCr = zeros(Height/2,Width/2);
case '4:2:2'
xqCb = zeros(Height,Width/2);
xqCr = zeros(Height,Width/2);
end
acBitsCb = 0;
dcBitsCb = 0;
acBitsCr = 0;
dcBitsCr = 0;
end
% Compute the bits for the Y component
for m = 1:N:Height
for n = 1:N:Width
t = y(m:m+N-1,n:n+N-1) - 128;
Yy = dct2(t); % N x N 2D DCT of input image
% quantize the DCT coefficients
temp = floor(Yy./(Qscale*QstepsY) + 0.5);
% Calculate bits for the DC difference
if n==1
DC = temp(1,1);
dcBitsY = dcBitsY + jpgDCbits(DC,'Y');
else
DC = temp(1,1) - DC;
dcBitsY = dcBitsY + jpgDCbits(DC,'Y');
DC = temp(1,1);
end
% Calculate the bits for the AC coefficients
ACblkBits = jpgACbits(temp,'Y');
acBitsY = acBitsY + ACblkBits;
% dequantize & IDCT the DCT coefficients
xqY(m:m+N-1,n:n+N-1)= idct2(temp .* (Qscale*QstepsY))+ 128;
end
end
% If the input image is a color image,
% calculate the bits for the chroma components
if Depth >1
if strcmpi(SamplingFormat,'4:2:0')
EndRow = Height/2;
else
EndRow = Height;
end
for m = 1:N:EndRow
for n = 1:N:Width/2
t1 = Cb(m:m+N-1,n:n+N-1) - 128;
t2 = Cr(m:m+N-1,n:n+N-1) - 128;
Ycb = dct2(t1); % NxN2DDCTof Cbimage
Ycr = dct2(t2);
temp1 = floor(Ycb./(Qscale*QstepsC) + 0.5);
temp2 = floor(Ycr./(Qscale*QstepsC) + 0.5);
if n==1
DC1 = temp1(1,1);
DC2 = temp2(1,1);
dcBitsCb = dcBitsCb + jpgDCbits(DC1,'C');
dcBitsCr = dcBitsCr + jpgDCbits(DC2,'C');
else
DC1 = temp1(1,1) - DC1;
DC2 = temp2(1,1) - DC2;
dcBitsCb = dcBitsCb + jpgDCbits(DC1,'C');
dcBitsCr = dcBitsCr + jpgDCbits(DC2,'C');
DC1 = temp1(1,1);
DC2 = temp2(1,1);
end
ACblkBits1 = jpgACbits(temp1,'C');
ACblkBits2 = jpgACbits(temp2,'C');
acBitsCb = acBitsCb + ACblkBits1;
acBitsCr = acBitsCr + ACblkBits2;
% dequantize and IDCT the coefficients
xqCb(m:m+N-1,n:n+N-1)= idct2(temp1 .* (Qscale*QstepsC))+ 128;
xqCr(m:m+N-1,n:n+N-1)= idct2(temp2 .* (Qscale*QstepsC))+ 128;
end
end
end
%
mse = std2(y-xqY)
snr = 20*log10(std2(y)/mse);
psnr = 10*log10(255*255/mse);

if Depth == 1
TotalBits = acBitsY + dcBitsY;
% title(['JPG compressed ' '@ ' num2str(TotalBits/(Height*Width)) ' bpp']);
axes(handles.jpeg);

imshow(xqY,[]);
disp('Thong so cua JPEG')
sprintf('SNR = %4.2f\n',snr)
sprintf('PSNR= %4.2f\n',psnr)

str = sprintf('Bit rate = %4.2f bpp\n\nMSE = %4.2f\n\nSNR = %4.2fdB\n\nPSNR= %4.2fdB\n\n',TotalBits/(Height*Width),mse,snr,psnr);
set(handles.jpegInfo,'String',str);

else
TotalBits = acBitsY + dcBitsY + dcBitsCb +...
acBitsCb + dcBitsCr + acBitsCr;
c1 = imresize(xqCb,[Height Width],'cubic');
c2 = imresize(xqCr,[Height Width],'cubic');
mseb = std2(A(:,:,2)-c1)%cai nay cua Cb
snrb = 20*log10(std2(A(:,:,2))/mseb)
psnrb =10*log10(255*255/mseb)
msec = std2(A(:,:,3)-c2)%cai nay cua Cr
snrc = 20*log10(std2(A(:,:,3))/msec)
psnrc = 10*log10(255*255/msec)

sprintf('PSNR(Cb) = %4.2fdB\tPSNR(Cr) = %4.2fdB\n',psnrb,psnrc)

sprintf('SNR(Cb) = %4.2fdB\tSNR(Cr) = %4.2fdB\n',snrb,snrc)
xq(:,:,1) = xqY;
xq(:,:,2) = c1;
xq(:,:,3) = c2;
%title(['JPG compressed ' '@ ' num2str(TotalBits/(Height*Width)) ' bpp'])
axes(handles.jpeg);
imshow(ycbcr2rgb(uint8(round(xq))));

str = sprintf('Bit rate = %4.2f bpp\n\nMSE(Y) = %4.2f\nMSE(Cr) = %4.2f\nMSE(Cb) = %4.2f\n\nSNR(Y) = %4.2fdB\nSNR(Cb) = %4.2fdB\t\nSNR(Cr) = %4.2fdB\n\nPSNR(Y)= %4.2f\nPSNR(Cb) = %4.2fdB\t\nPSNR(Cr) = %4.2fdB\n',TotalBits/(Height*Width),mse,mseb,msec,snr,snrb,snrc,psnr,psnrb,psnrc);
set(handles.jpegInfo,'String',str);
end

